package com.akasiyanik.bsu.coursework.methods;

import static com.akasiyanik.bsu.coursework.utils.MatrixUtils.*;

import com.akasiyanik.bsu.coursework.equations.SteadyingEquation;
import com.akasiyanik.bsu.coursework.methods.rungekutta.AuxRungeKuttaMethod;
import com.akasiyanik.bsu.coursework.utils.MatrixUtils;

/**
 * @author: akasiyanik
 */
public class SteadyingProcess {

    private static double EPS;
    private SteadyingEquation steadyingEquation;
    private double w;
    private AuxRungeKuttaMethod auxMethod;

    public SteadyingProcess(AuxRungeKuttaMethod auxMethod, double EPS, double w, SteadyingEquation steadyingEquation) {
        this.auxMethod = auxMethod;
        this.EPS = EPS;
        this.w = w;
        this.steadyingEquation = steadyingEquation;
    }

    public double[] getY() {
        int s = steadyingEquation.getS();
        int n = steadyingEquation.getN();
        double[] y0 = steadyingEquation.getY0();
        int dim = s * n;

        //составляем начальный вектор Y, состоящий из s векторов y0, размерности n
        double[] Y = new double[dim];
        for (int i = 0; i < s; i++) {
            for (int j = 0; j < n; j++) {
                Y[i * n + j] = y0[j];
            }
        }
        double err;
        int iterationNumber = 0;
        do {
            double[] Y1;
            if (iterationNumber % 10 == 0) {
                 Y1 = doSteadyingIterationWithClarifying(Y);
            } else {
                Y1 = doSimpleSteadyingIteration(Y);
            }
            err = calculateErr(Y, Y1);
//            System.out.println("POST In method Err = " + err);
            Y = Y1;
            iterationNumber++;
        } while (err >= EPS);

        return Y;
    }

    private double[] doSteadyingIterationWithClarifying(double[] Y) {
        double[] Y1 = F(Y);
//        System.out.println("---------------- Y1 ----------------");
//        printArray(Y1);
//        double err = calculateErr(Y, Y1);
//        System.out.println("\nPRE In method Err = " + err);
        if (checkingConvergence()) {
            return clarify(Y, Y1);
        } else {
            return Y1;
        }
    }

    private double[] doSimpleSteadyingIteration(double[] Y) {
       return F(Y);
    }

    private boolean checkingConvergence() { // проверка сходимости
        // todo akasiyanik
        return true;
    }

    private double[] clarify(double[] Y, double[] Y1) {
        double[] r_k = steadyingEquation.r(Y);
        double[] r_k1 = steadyingEquation.r(Y1);
        double lambda_m = getEigenvalueApproximation(r_k, r_k1);
        Y1 = MatrixUtils.subtract(Y1, MatrixUtils.scalarMultipy(w / lambda_m, r_k1));
        return Y1;
    }

    //среднее значение центрального квантиля

    private double getEigenvalueApproximation(double[] r_k, double[] r_k1) { //getting lambda_m
        double derRin0 = auxMethod.getDerivativeRin0();
        double[] v = MatrixUtils.perComponentDivisionWithZerosInFraction(r_k1, r_k);
        //todo akasiyanik instead of average it is required to use average of central quantile
        double t = MatrixUtils.getMedianWithoutZeros(v);
//        System.out.println("!!!median: " + t);
        double eigenvalueApproximation = (t - 1.0) / derRin0;
//        System.out.println("!!!eigenvalueApproximation: " + eigenvalueApproximation);

        return eigenvalueApproximation;
    }

    private int getMostHeavyComponentNumber(double[] Y, double[] Y1) {
        double maxDiff = 0.0;
        int mostHeavyComponent = 0;

        for (int i = 0; i < Y.length; i++) {
            if (Math.abs(Y1[i] - Y[i]) >= maxDiff) {
                maxDiff = Math.abs(Y1[i] - Y[i]);
                mostHeavyComponent = i;
            }
        }
        return mostHeavyComponent;
    }


    private double calculateErr(double[] arr1, double[] arr2) {
        double max_err = 0.0;
        for (int i = 0; i < arr1.length; i++) {
            double err = Math.abs(arr1[i] - arr2[i]);
            if (err > max_err) {
                max_err = err;
            }
        }
        return max_err;
    }

    private double[] F(double[] Y0) {
        int n = auxMethod.getN();
        double[][] betta = auxMethod.getBetta();
        double[] Y = Y0;
        double gamma = auxMethod.getGamma();

        if (auxMethod.getSigma() % 2 == 0) { //четное сигма
            for (int m = 0; m < n; m++) {
                double[] sum1 = scalarMultipy(betta[m][0], steadyingEquation.r(Y));
                double[] wr = scalarMultipy(w, steadyingEquation.r(Y));
                double[] sum2 = scalarMultipy(betta[m][1], steadyingEquation.r(add(Y, wr)));
                Y = add(Y, scalarMultipy(w, add(sum1, sum2)));
            }
        } else { // нечетное сигма
            for (int m = 0; m < n - 1; m++) {  // n-1 двухстадийные методы
                double[] sum1 = scalarMultipy(betta[m][0], steadyingEquation.r(Y));
                double[] wr = scalarMultipy(w, steadyingEquation.r(Y));
                double[] sum2 = scalarMultipy(betta[m][1], steadyingEquation.r(add(Y, wr)));
                Y = add(Y, scalarMultipy(w, add(sum1, sum2)));
            }
            Y = add(Y, scalarMultipy(w * gamma, steadyingEquation.r(Y)));   //n-ый метод - метод Эйлера
        }
        return Y;
    }

//    private double[][] calculateK(double[] Y) {
//        int sigma = auxMethod.getS();
//        int n = steadyingEquation.getN();
//        int s = steadyingEquation.getS();
//        double[][] a = auxMethod.getA();
//        double[][] K = new double[sigma][s * n];
//        for (int p = 0; p < sigma; p++) {
//            double[] sum = new double[s * n];
//            for (int q = 0; q <= p - 1; q++) {
//                sum = add(sum, scalarMultipy(a[p][q], K[q]));
//            }
//            K[p] = steadyingEquation.r(add(Y, scalarMultipy(w, sum)));
//        }
//        return K;
//    }


}
